[1]:

import pylab as pp
import numpy as np
import pandas as pd

from scipy import integrate, interpolate, optimize
from scipy.integrate import odeint

from bokeh.palettes import brewer
from bokeh.plotting import figure, show
from bokeh.io       import output_notebook

output_notebook()
Loading BokehJS ...

Dados artificiais

Uma prática interessante na análise de dados é testarmos se o algoritmo proposto e principalmente sua estrutura está consistente. Podemos então utilizar dados artificiais para avaliar se o processo de ajuste está estruturado corretamente.

[2]:


# Lendo o arquivo de dados no formato 'filename.csv'
data = pd.read_csv("./PG_IMT/DadosEpidemia/SIRpd")
# Preview das cinco primeiras linhas
data.head()


[2]:
S I R
0 499.000000 1.000000 0.000000
1 498.949042 1.040746 0.010213
2 498.896013 1.083146 0.020841
3 498.840830 1.127267 0.031903
4 498.783405 1.173179 0.043415 
[3]:


s_array = data[["S", "I", "R"]].to_numpy()

Sd = s_array[:,0]
Id = s_array[:,1]
Rd = s_array[:,2]

Gerando ruído gaussiano

[4]:


Sdn = np.random.normal(0, np.mean(Sd)/10, len(Sd))
Idn = np.random.normal(0, np.mean(Id)/5,  len(Sd))
Rdn = np.random.normal(0, np.mean(Rd)/10, len(Sd))

Sd = Sd + Sdn
Id = Id + Idn
Rd = Rd + Rdn

Visualizando os dados

[5]:


# Visualizando a evolução da Epidemia - S(t), I(t) e R(t)
t = np.linspace(0,len(Sd),len(Sd))

TOOLS="hover,crosshair,pan,wheel_zoom,zoom_in,zoom_out,box_zoom,undo,redo,reset,tap,save,box_select,poly_select,lasso_select,"

p = figure(tools=TOOLS, plot_width=600, plot_height=400)

p.scatter(t, Sd, legend_label="Suscetíveis - dados",
          radius=3.8, fill_color="#ffd885", fill_alpha=0.6, muted_color="#ffd885", muted_alpha=0.2, line_color=None)
p.scatter(t, Id, legend_label="Infectados - dados",
          radius=3.8, fill_color="#de425b", fill_alpha=0.6, muted_color="#de425b", muted_alpha=0.2, line_color=None)
p.scatter(t, Rd, legend_label="Removidos - dados",
          radius=3.8, fill_color="#99d594", fill_alpha=0.6, muted_color="#99d594", muted_alpha=0.2, line_color=None)

p.grid.grid_line_alpha = 0
p.ygrid.band_fill_color = "olive"
p.ygrid.band_fill_alpha = 0.1
p.yaxis.axis_label = "Indivíduos"
p.xaxis.axis_label = "Dias"
p.legend.click_policy = "mute"
p.legend.items.reverse()

show(p)

O problema

O conjunto de equações diferenciais que caracteriza o modelo é descrito abaixo. No modelo $:nbsphinx-math:beta `- :nbsphinx-math:text{representa a taxa de transmissão ou taxa efetiva de contato}` $ e \(r - \text{a taxa de remoção ou recuperação.}\)

\[\begin{split}\begin{split} \frac{dS(t)}{dt} & = -\beta S(t) I(t) \\ \frac{dI(t)}{dt} & = \beta S(t) I(t) - rI(t) \\ \frac{dR(t)}{dt} & = r I(t) \end{split}\end{split}\]

Gostaríamos de identificar quais parâmetros \(\beta\) e \(r\) resultam num melhor ajuste do modelo para os dados de S,I e R

[6]:


def SIRmodel(y, t, Beta,r):
    S, I, R = y
    Sdot = -(Beta * S * I)
    Idot = (Beta * S * I)  - r * I
    Rdot = r * I
    return Sdot, Idot, Rdot

# Resolução da simulação - Escala temporal (dias)

Obtendo \(y_s(\theta,k) = [S \; I \: R]\)

O trecho a seguir retorna os valores sintetizados \(y_s(\theta,k) = [S \; I \: R]\) representa o dado sintetizado a partir de um modelo sintetizado para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\). A partir de uma condição inicial \(y_0\).

[7]:


def SIRsim(y0,t,theta):
    Beta = theta[0]
    r = theta[1]
    ret = integrate.odeint(SIRmodel,y0,t,args=(Beta,r))
    S, I, R = ret.T
    return S, I, R

Condições inicias - \(y_0\) e \(\theta_0\)

[8]:


# Tamanho da populção - N
N = 500

# Valores iniciais
I0, R0 = 1, 0
S0 = N - I0

# Vetor de condições iniciais

y0 = S0, I0, R0

# Beta -  taxa de contato,
# r - taxa média de recuperação (in 1/dia).

theta0 = [1e-4,1e-2] # valores iniciais

# Definição do conjunto de equações diferencias não lineares que formam o modelo.

t = np.linspace(0, 1000, 1000)

Estimação de parâmetros

Para estimarmos os parâmetros do modelo \(\mathbf{\beta}\) e \(\mathbf{r}\), vamos utilizar inicialmente o método de mínimos quadrados. Podemos então formular o problema a partir da Equação abaixo. Na Equação \(y_m(k)\) representa o dado real em cada amostra \(k\); \(y_s(\theta,k)\) representa o valor estimado a partir da simulação do modelo para uma determinada amostra \(k\) e \(\theta\) representa o vetor ed parâmetros \(\theta = [ \beta \; \; r]^T\).

\[min_{\theta}= \sum_{k=1}^{K}(y_m(k) - y_s(\theta,k))^2\]

A equação formula a pergunta: quais os valores de \(beta\) e \(r\) que minizam o erro quadrático quando comparados com os dados reais.

[9]:


def ErroQuadratico(Sd,Id,Rd,y0,t,theta0):
    """ function to pass to optimize.leastsq
        The routine will square and sum the values returned by
        this function"""
    [S,I,R] = SIRsim(y0,t,theta0)
    erroS = S - Sd
    erroI = I - Id
    erroR = R - Rd
    EQ = np.concatenate([erroI,erroR])
    return EQ

def objetivo(p):
    return ErroQuadratico(Sd,Id,Rd,y0,t,p)

Minimização da função custo

[10]:


(c,kvg) = optimize.leastsq(objetivo,theta0)
print(c)


[9.93827391e-05 9.96628326e-03]

Visualização

[11]:


[Sa,Ia,Ra] = SIRsim(y0,t,c)

# Visualizando a evolução da Epidemia - S(t), I(t) e R(t)
p.line(t, Sa, legend_label="Suscetíveis", color="#f9bd3d", line_width=3)
p.line(t, Ia, legend_label="Infectados",  color="#f4193c", line_width=3)
p.line(t, Ra, legend_label="Removidos",   color="#45c83a", line_width=3)

show(p)